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Abstract 

We study the emergence of magnetic flux from the near-surface layers of the solar convection zone 
into the photosphere. To model magnetic flux emergence, we carried out a set of numerical radiative 
magnetohydrodynamics simulations. Our simulations take into account the effects of compressibility, energy 
exchange via radiative transfer, and partial ionization in the equation of state. All these physical ingredients 
are essential for a proper treatment of the problem. Furthermore, the inclusion of radiative transfer allows 
us to directly compare the simulation results with actual observations of emerging flux. We find that the 
interaction between the magnetic flux tube and the external flow field has an important influence on the 
emergent morphology of the magnetic field. Depending on the initial properties of the flux tube (e.g. field 
strength, twist, entropy etc.), the emergence process can also modify the local granulation pattern. The 
emergence of magnetic flux tubes with a flux of 10 19 Mx disturbs the granulation and leads to the transient 
appearance of a dark lane, which is coincident with upflowing material. These results are consistent with 
observed properties of emerging magnetic flux. 

1 Introduction 

Solar magnetic fields on the surface of the Sun exist and evolve over a wide range of length- and time- 
scales. The most prominent magnetic features on the solar surface are sunspots. In addition, there is a 
whole hierarchy of magnetic features including pores, micropores, plages and faculae. An active region is 
an extended bipolar c onfiguration o n the solar surface resulting from the emergence of magnetic fields from 
the convection zone (|Parkerl Il955h . In term s of the amount of m agnetic flux in each polarity, there is 
a continuous spectrum of active region sizes l|Hagenaar et all I2003T I . Although a partitioning of the flux 
spectrum for the sake of classification may seem somewhat arbitrary, it allows us to conveniently refer to 
active regions of different sizes. Large active regions have polarities containing a flux of > 5 x 10 21 Mx 
and contain sunspots. In large active regions, the magnetic flux is shared amongst the whole spectrum of 
magnetic features. Large active regions have lifetimes of up to months. Small active regions, which contain 
a flux of 1 x 10 20 — 5 x 10 21 Mx in each polarity, may consist of pores and smaller magnetic features but 
lack sunspots. Small active regions may persist for up to days to weeks. Ep hemera l activ e regions have even 
less flux (3 x 10 18 — 1 x 10 20 Mx), and have lifetimes of only hours to days l)Zwaanl . ll98^ . Often, ephemeral 
active regions are simply referred to as ephemeral regions. 

The characteristic timescale for the emergence of flux is typically only a small fraction of the lifetime 
of an active region. For instance, almost all the flux of a large active region emerges within the first ~ 4 
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days of their development l|Zwaanl . 1 19851 : lHagenaaij , l200ll ). Observations of the birth of an active region 
(AR) indicate that the total flux in each polarity of the AR is not the consequence of the emergence of a 
coherent, monolithic flux bundle. Rather, it builds up as the result of many small flux bundles emerging 
simultaneously or in succession. A n Emerging Flux Region (EFR) is the area on the solar surface where 
these emergence events take place ijZirinl [l972l : IZwaanl . Il985h . The onset of the birth of an active region is 
characterized by the appearance of a compact and very bright plage. The plage consists of magnetic flux 
elements of opposite polarity, which move apart at an initial velocity of about 2 km s _1 . New flux continues 
to emerge near the polarity inversion line. The orientation of the emerging field is not random. Rather, it is 
roughly aligned along the axis connecting the two polarities. As a result of emergence, flux is accumulated in 
both polarities. If sufficient flux has emerged, pores and possibly sunspot s appe ar. These tend to be formed 
near the leading and following edges of the expanding plage (see IZwaanl . 1 19851 . and references therein) . 

The granulation pattern in an EFR may appear different than that of the quiet Sun. In quiet-Sun granu- 
lation, bright granules correspond to upflows whereas the dark intergranular boundaries consist of downflow 
lanes and vertices. Transient dark alignments in the central part of an EFR ha ve been detected in continuum 
images and also in maps of the core intensity of photospheric spectra l lines (jBrav. Louehhead. fc DurrantL 
11984 ; iBrants fe Steenbeekl Il985l : IZwaan . 1985; S trous fc Zwaanl . 1 19991 ). The darkenings are roughly aligned 
along the axis connecting the two polarities of the active region and typically last about 10 minutes. In the 
continuum, they are darker than the intergranular boundaries but the spectral lines show Doppler shifts cor- 
responding to upward velocities exceeding 0.5 km s _1 . By estimating the d iameter and strength of eme rging 
flux bundles that lead to the appearance of dark alignments in an EFR, IBrants fc Steenbeekl (| 1985h esti- 
mated that each bundl e conta ins a longitudinal flux of about 10 19 Mx. This value is consistent with the 
estimate given bv lBornl (|l9T4h . who measured the total flux of an active region and counted the number of 
arch filaments in Ha. 

On the basis of the observations previously summarized, IZwaanl |l97Sl . Il985l ) constructed a heuristic 
model of an EFR at and below the solar surface. In his scenario, flux emerges as a collection of arched 
flux tubes rising through the convection zone. Deep in the convection zone, the flux tubes connect to an 
underlying coherent large flux tube. Near the surface, however, they are separated from each other. In his 
model, the transient dark alignments in the observations are locations of the apices of arched flux tubes 
emerging at the surface, which appear dark because the horizontal magnetic field suppresses turbulent heat 
exchange. After the horizontal top of a tube has emerged, the photospheric footpoints of the tubes separate 
and the field at the footpoints become increasingly vertical. The coalescence of the vertical flux elements in 
each polarity may lead to the formation of pores and sunspots if sufficient flux is accumulated. 

Su bsequent observations of EFRs are consistent with various aspects of Zwaan's heuristic model. Istrous fc Zwaanl 
1 19991 ) performed a statistical analysis of over two hundred emergence events in a single EFR. They found 
that emerging flux is characterized by the transient appearance of dark alignments between the polarities. 
Often, they observed the appearance of faculae at the ends of th e dar k alignments . The faculae are typi- 
cally associated with downflows. More recently, iLites et all (|l998h and lKubo et~aH (|2003T ) have carried out 
observations of the full Stokes vector in emerging flux regions and report the detection of horizontal fields 
with strengths of about 200 — 600 G. As the newly emerged flux moves away from the emergence site, the 
fields become vertical and obtain strengths in excess of 1000 G. 

The aim of this study is to understand the process of magnetic flux emergence from the convection zone 
into the photosphere by means of numerical simulations. At present, it is not computationally feasible to 
carry out numerical simulations that include all the relevant physics for the emergence of magnetic flux 
through all layers of the solar convection zone and atmosphere. To progress, different aspects of the problem 
have to be addressed separately. 

One line of previous research has focused on studying the rise of buoyant magnetic flux tu bes in a stratified 
layer by means of 2D and 3D MHD simulations. Most of these simulations (see, for example. ISchusslerl . 
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Cheung et all l2006h start from a buoyant magnetic flux tube embedded in an initially static atmosphere 



Porch et all (|200ll ) was the first to study the interaction of magnetic flux tubes with convective flows. They 



found that the convective motion of the flow external to the tube en hances the loss rate of flux from the 
tube. IFan. Abbett. fc Fisher] (120031 ) and lAbbett. Fisher, fc Fan! ()2000h have also studied the interaction of 
magnetic flux tubes and convective flow. They found that flux tubes with field strengths exceeding the 
equipartition value with respect to the kinetic energy density of the convective motions are able to resist 
deformation by the convective flow, while significantly weaker flux tubes are severely distorted. 

Another line of work has focused on studying the dynamics of toroidal magnetic flux tubes rising in 
the solar conve ction zone. Numerical simulat i ons o f rising toroidal flux tubes under the thin flux tube 
approximation (R oberts fc Webbl Il978l ; ISpruitl Il98ll ) have been successful at reproducing the large-scale 
properties of active regions, such as thei r emergence latitudes, their tilt angles as we l l as the asymmetries 
between lea ding and following polarities l|D'Silva fc Choudhuril [l993l ; IFan et all [l993. 19941; ICaligari et all . 
1 19951 . 1 19981 ). When the apices of magnetic loops reach a depth of about 10 — 20 Mm below the photosphere, 
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their cross-sections have expanded so much that the thin flux tube approximation no longer applies. This 
limitation means that fully multi-dimensional MHD simulations must be carried out in order to study the 

details of the emergence process. 

Beginning with the work of iForbes fc Priestl (|l984h and IShibata et all l|l989l) . a variety of flux emer- 
gence simulations have been carried out. The focus of most studies in the literature has been on the 
evolution of the emerging magnetic field in the corona (fFanL l200ll: iMaeard. I200I; iMaeara fc Longcopel. l2003| : 
Fan fc Gibsonl.l2004lAbbett fc Fisherl,l2003] ; lArchontis et al.l.l2004l2005l;lTorok fc Klieml,l2005l ; llsobe et"aT 



2005l;lGalsgaard et al.l . l2005l : lLeake fc Arber[|2006l : lMurrav et al.l . l2006l : iGalsgaard et alll2007h . Recently. iMagaral 



(2006) carried out MHD simulations of flux emergence to study the evolution of the coronal and photospheric 
field of an emerging bipolar region. Like the aforementioned references, his work is relevant for the large 
scale behaviour of active regions and ignores the effects of convection as well as radiative energy exchange. 

In this study, we examine in detail the effects of convective flows and radiative transfer on emerging 
magnetic flux. To this end, we restrict our attention to flux emergence from the near-surface part of the 
convection zone into the photosphere. The article is structured as follows. In Sect. [2] we present the system 
of radiative MHD equations, the numerical methods used to solve this system as well as the setup of the 
simulations. In Sect. [3] we discuss the dynamics of rising flux tubes with special emphasis on their interaction 
with the convective motions. In Sect. [4] we discuss the observational signatures of emerging flux as obtained 
from our simulations and compare them with observational studies in the literature. Finally, we close the 
paper with a discussion in Sect. 



2 Governing equations, numerical methods and simulation 
setup 

2.1 Governing equations 

The system of radiative MHD equations (in Gaussian cgs units) consists of 
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where g is the fluid mass density, v the fluid velocity, B the magnetic field, p to t = p + \B\ 2 /Stv the sum of 
the gas and magnetic pressures, and e = ge + ^p\v\ 2 + |-B| 2 /87r the sum of the internal, kinetic and magnetic 
energy densities per unit volume. The symbol ® denotes a tensor product and 1 is the 3x3 identity matrix. 
T is the gas temperature and K the thermal conductivity, r represents the viscous stress tensor. Q ra d is 
the radiative heating rate and g = —2.74 x 10 4 i cm s -2 is the gravitational acceleration. 

Equations ([2]) and ((3| express, respectively, the principles of mass conservation, momentum balance 
and energy conservation, while Eq. @ is the induction equation, which describes the evolution of a magnetic 
field in the plasma. In this equation, 77 = c 2 /Aira is the magnetic diffusivity, c the speed of light and a the 
electrical conductivity. 

The MURaM code was used to carry out the radiative MHD simulations in this study. The simulations 
presented here used grey radiative transfer (i.e. 1 frequ e ncy b in). For details about th e numerical methods 
implemented in MURaM, we refer the reader to lVogleri l|2003h and lVogler et all (|2005h . 



2.2 Boundary conditions 

Periodic boundary conditions are imposed at t he vertical boundar ies. The open boundary condition at the 
lower boundary is identical to the one used bv lVogler et~ai1 (|2005l ). 

The top boundary condition has been modified to allow for the advection of magnetic field through 
the top boundary. This is important because we do not want the emerged magnetic field to be artificially 
trapped in the photosphere. In principle, the following conditions 

dv x dv y dv z 
oz oz az 
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Label 


B (G) 
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$ (10 19 Mx) 


*SH ( 1Qiy Mx ) 


Emergence fraction (%) 


Ul 


8500 


0.5 


1 


7.7 


67 


U2 


8500 


0.25 


1 


4.7 


08 


U3 


8500 


0.1 


1 


3.1 


44 


U4 


8500 





1 


2.6 


30 


U5 


8500 


0.75 


1 


11 


62 


U6 


2500 


0.5 


0.31 


0.9 


12 


U7 


5250 


0.5 


0.66 


3.0 


39 


U8 


7000 


0.5 


0.86 


5.6 


49 



Table 1: Initial properties of the individual magnetic flux tubes and emerged flux of the simulation runs. The 
columns give (from left to right) the label of the run, the magnitude of the longitudinal field at the tube axis 
(Bo), the twist parameter (A), the total longitudinal flux ($0)1 the maximum unsigned magnetic flux ($J^c) 
crossing the z — plane (i.e. base of the photosphere) within a time period of 20 min after the start of 
the simulation, and the percentage of tracer fluid elements in the initial flux tube that manage to rise to the 
photosphere (z — 0) within the same time period. In all cases, the characteristic radius of the initial flux tube 
is R = 200 km. 



should suffice for a smooth outflow. This condition, however, is independent of whether the mass flux through 
the upper boundary is appropriate or unrealistically high. In ord er to keep the ma s s flux at appropriate 
levels, additional constraints must be applied. To this end, we follow l&tein fc Nordlundl (|l998l ) and implement 
a fiducial layer above the top boundary of our simulation domain. 

The magnetic field above the upper boundary of the domain is matched to a potential field. This requires 
an extrapolation of the magnetic field in the uppe r most domain layer into the ghost cells at each time step. 
For further details, we refer the reader to ICheuna (|2006l ). 

2.3 Initial conditions 

We have carried out a number of simulation runs, each modelling the emergence of an individual flux tube 
initially embedded in the near-surface layers of the convection zone. In order to prepare the initial condition 
for the flux emergence simulations, a 3D dynamic model of the near-surface layers of the convection zone and 
photosphere is required. We obtained such a model by means of a non- magnetic simulation with the MURaM 
code. The details of the corresponding simulation setup are given in ICheung. Schiissler. fc Moreno-Insertisl 
(|2007h . 

The simulation domain of our model atmosphere has dimensions 24 x 12 x 2.3 Mm 3 . The level z — 
corresponds to the mean geometrical height where the continuum optical depth at 5000 A is unity (i.e. 
T5000 = 1) and is located at a height of 1.85 Mm above the bottom boundary of the domain. We implant 
a horizontal, twisted magnetic flux tube with the axis of the tube located at [y,z] = [6, —1.35] Mm into a 
snapshot of the 3D hydrodynamic simulation. In all subsequent discussion, t = refers to the moment when 
the flux tube is introduced. The longitudinal and transverse (azimuthal) components of the magnetic field 
have the form 

B,(r) = B exp(-r 2 /^o), (6) 

Be(r) = (7) 
Ho 

where r £ [0, 2R$] is the radial distance from the tube axis. The longitudinal flux carried by the tube is 
$0 = / BidS ~ 0.987ri?g_B , so that R can be viewed as a characteristic 'tube radius'. The dimensionless 
twist parameter A specifies the relative strength of the transverse component of the magnetic field with 
respect to the longitudinal component. 

Let Vorig be the original velocity distribution in the non-magnetic domain. In order to start with an 
almost stationary flux tube without excessively perturbing the convective flow, we impose the following 
initial velocity distribution within the tube: v tu bc (r) = (l-e~ r /R o)v orig . With increasing radial distance 
from the axis, the fluid velocity approaches the original velocity distribution of the external convecting flow. 

For a given magnetic field distribution in the tube, the internal gas pressure distribution is specified by 
requiring that the divergence of the combined stress tensor (gas pressure + Maxwell) be identical to that 
in the non-magnetic domain (the divergence of the viscous stress tensor is not identical because we have 
modified the velocity distribution in the magnetic flux tube). Having done so, one is still free to choose the 
distribution in the tube of one of the following thermodynamic properties: mass density, specific entropy, 
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or temperature. In all of the simulations presented in this paper, the initial flux tube was imparted with a 
uniform specific entropy distribution (r < 2Ro). The value of the specific entropy was chosen as its average 
value in the convective upflows at the initial depth of the tube. 



3 Dynamics of subsurface rise and appearance at the photo- 
sphere 

3.1 The rise of a buoyant magnetic flux tube in the presence of convective 
flows 

Convection plays an important role in determining the properties of emerging flux. Near-surface con- 
vectio n in the Sun consists o f disjointed upwellings separated by a network of downflow lanes and ver- 
tices (jStein fc Nordluncfl . Il989f ). An initially horizontally flux tube embedded in this setting will encounter 
both upflows and downflows. Whereas the upwellings aid the emergence of some segments of the tube, the 
downflows impede the rise of other parts. Given such an initial configuration, under what circumstances can 
we expect the tube to emerge coherently without severe distortion by the downflows? In other words, what 
does it mean to have a 'strong' flux tube? 

We expect that the flux tube can resist severe distortion by the flow if the field strength satisfies Bo > -B cq , 
where B cq is the field strength at which the magnetic energy density and the kinetic energy density of the 
external flow are in equipartition. The value of the equipartition field strength depends on depth. At the 
original location of the flux tube (z — —1.35 Mm), the ambient density is g = 4.2 x ICP 6 g cm~ 3 . The 
typical vertical velocity of downflowing material at this depth is 4 — 8 km s _1 . For this range of velocities, 
the corresponding range of equipartition field strengths is B eq = 2900 — 5800 G. Taking the typical values 
at the surface (g = 2.6 x 10 -7 g cm -3 , v — 2 — 4 km s _1 ), the equipartition field strength there has much 
lower values of B cq = 450 - 700 G. 

Consideration of the force balance on the tube allows us to determine a related, but more stringent 
criterion. Suppose that the flux tube has uniform magnetic field strength, Bo, and uniform density deficit, 
Ag, with respect to the surroundings. The buoyancy force experience d by the flux tube overcomes the drag 
force exerted by a downflow of speed v given th at (see, for example. IParkeil Il979l ; [Moreno-Insertis! . Il983l : 
ISchiisslerl . Il98l ISchiissleil Il987l : iFan et all l2003h 

B.af^y'^V'V (8) 



where H p is the local pressure scale height, Cd the drag coefficient (of order unity) and 71 the first adiabatic 
exponent. The factor 2Cr>7i/7r is of order unity. The radius is important because the drag force (per unit 
length) on the flux tube is proportional to the tube radius while the buoyancy force (also per unit length) is 
proportional to the square of the tube radius and to the square of the field strength. This means that, for a 
given field strength, a thick tube tends to be dominated by buoyancy while a sufficiently thin tube is passively 
dragged by the surrounding flows. For example, we see from the criterion given by Eq. ([8} that a tube with 
a radius Rq — 0.25H P must have a field strength of at least 2B cq in order to have sufficient buoyancy to rise 
against the drag of the downflows. Since the downflows are localized, the magnetic curvature of a deformed 
flux tube also acts against the drag force. Like the buoyancy force, the curvature force is also proportional 
to the square of the magnetic field strength. 

Table I gives on overview of the various simulation runs carried out. We start the discussion of the results 
with the 'weak' flux tube considered in run U6. The initial field strength at the tube axis is Bo — 2500 G, 
corresponding to f3 — 8irp/B 2 ~ 22 there. The characteristic radius of the tube is Ro = 200 km ~ 0AH p . 
The tube is initially twisted, with A = 0.5. The longitudinal magnetic flux threading a cross-section of 
the tube is 3.1 x 10 18 Mx. Criterion (|8]) tells us that the tube must have at least a field strength of 
(H p /R ) 1/2 B etl = 4600 - 9200 G in order to rise against the downflows. Since the tube has only a central 
field strength of Bo = 2500 G (the average field strength over its cross-section is even smaller), we expect 
the evolution of the flux tube to be dominated by the convective flows. Fig. [1] shows a time sequence of 
3D isosurfaces of \B\ = 400 G. In all three panels, the blue (red) colour-coding indicates that the vertical 
velocity at the isosurface is upwards (downwards). Already at 2.2 min after the start of the simulation, the 
flux tube has been significantly distorted by the convective flows. In the subsequent snapshots, the upflows 
have advected segments of the tube upwards, whereas the downflows have pinned down other segments below 
the surface. The shape of the flux tube has developed to a form resembling a sea serpent. The behaviour of 
the flux tube in the presence of asymmetrical u p- and down-flows in compressible convection is reminiscent 
of the turbulent pumping mechanism studied bv lTobias etafl l|200ll ). 
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Figure 1: The passive evolution of a weak magnetic flux tube (run U6, Bo = 2500 G) with the convective flow. 
Shown is a sequence of isosurfaces of \B\ = 400 G, on which the blue (red) colour-coding indicates upflows 
(downflows). Segments of the tube embedded in upflows are able to rise and emerge, whereas segments aligned 
with downflows are kept submerged. 



Figure [5] gives another illustration of how the flux tube in run U6 evolves passively with the convecting 
flow. The figure shows a sequence of snapshots of a vertical cross-section of the simulation domain at x = 12 
Mm (perpendicular to the axis of the initial flux tube) . The grey scale indicates the temperature distribution 
while the colour-coding represents the absolute field strength, The arrows indicate the components of 
the velocity field in the yz plane at x = 12 Mm. The purple line running near z = shows the level of 
optical depth unity in the continuum at a wavelength of 5000 A (rgooo = 1)- At £ = min, this particular 
segment of the tube is aligned with a downflow. In the following snapshot at t = 4.4 min, we find that this 
segment of the tube has been displaced downwards. In the final two snapshots (t = 8.8 min and t — 13.3 
min), we witness the tube segment being severely deformed by the shear in the velocity field at the interface 
between upflows and downflows. While the central part of the tube segment is pushed downwards, other 
parts of the tube segment are carried upwards. 

Our simulation runs with higher initial field strengths confirm that the ability of the convective flows to 
dominate the evolution of the tube lessens with increasing initial field strength. In run Ul, the initial field 
strength on the axis of the tube is Bo = 8500 G, which is within the range of B aq for the original depth of 
the tube and corresponds to a (3 ~ 1. The total longitudinal magnetic flux carried by the tube is 10 19 Mx. 
Fig- EH shows a time sequence of isosurfaces of \B\ = 700 G for this run. Again, the colour-coding indicates 
upwards or downwards vertical velocity. In this case, the flux tube is not so weak that it is simply advected 
by the convective flow. Nor is it sufficiently strong to be unaffected by the flows. The behaviour of the tube 
is in an intermediate regime between these two extremes. On the one hand, the convective flow is able to 
deform the tube. On the other hand, the bulk of the tube finally does overcome the impeding downflows 
and emerges at the surface. 

Figure |3 helps emphasize the previous point. At t = min, the tube segment shown is aligned with a 
downflow. In run U6, this downflow advects the tube downwards. The situation in run Ul is very different. 
While the downflow is able to divert the tube slightly towards the right, it is unable to keep the tube beneath 
the surface. Instead, the tube segment is sufficiently buoyant to overcome the drag by the downflow and 
finally emerge. 

3.2 Dependence of emerged flux on initial field strength and twist 
3.2.1 Unsigned flux at the photosphere 

The 'fraction' of the original flux tube which manages to emerge at the photosphere can be quantified in 
different ways. One measure is the maximum value reached in a given time interval by the total unsigned 
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Time = 0.0 min VB 




Time = 8.8 min 




3 4 5 6 7 8 9 

y (Mm) 300 



Figure 2: Time sequence of a vertical cross-section of the simulation domain at x = 12 Mm (perpendicular to the 
initial tube) for run U6. The tube segment is passively advected by the convecting flow. The greyscale indicates 
the temperature distribution and the color-coding indicates the distribution of the absolute field strength, \B\. 
The arrows indicate the projection of velocity vectors on the vertical cross-section. The purple line shows the 
level T5000 = 1- 
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Figure 3: Same as Fig. Q]but for a flux tube with higher initial field strength (run Ul, B = 8500 G). The 
isosurfaces here correspond to \B\ — 700 G. In this case, the convective flows do not completely control the 
dynamics of the flux tube. 



flux through the plane z = For a time interval of 20 min, this quantity is indicated in the fifth 

column of Table [T] Since a flux tube emerges as an undulated structure, the emergence actually consists of 
a collection of small-scale emergence events, each with a typical length-scale of about 1 Mm. This means 
that the unsigned vertical flux at the photosphere can be several times that of the initial longitudinal flux 
of the tube (<I?o). We identify two clear trends. Firstly, increases with increasing Bo for constant A. 

This simply reflects the fact that a magnetic flux tube with higher Bo contains more flux (for the same tube 
radius) and is more buoyant. Secondly, $j^x increases with A. This is due to the fact that the component 
of the magnetic field transverse to the original tube axis contributes to the unsigned flux at z = 0. 



3.2.2 Emergence of tracer fluid elements 

Another measure of the fraction of the tube that emerges can be obtained by considering selected (La- 
grangian) fluid elements as tracers of the evolution. To this end, we followed the trajectories of 250 such 
tracer fluid elements (referred to simply as 'tracers' in what follows). The first 50 of them were located 
equidistantly along the axis of the initial flux tube. The remaining tracers were selected similarly to lie 
along four lines parallel to the tube axis, each line at a distance of 50 km from the axis (one line lying above 
the axis, one lying below, and the other two lines to the left and right of the axis). We determined the 
fraction of tracers that manage to reach the photosphere [z > 0) within the time period < t < 20 min. 
The numerical values of the results indicated in the last column of Table [1] and plotted in Fig [5] For flux 
tubes with the same initial twist A, the percentage systematically increases with the initial field strength, 
Bo- As a reference case, we performed a non-magnetic run (i.e. Bo = 0) and found that 7% of the fluid 
elements emerged within 20 minutes. In the case of a weak magnetic flux tube (run U6, So = 2500 G), 
12% of the tracers emerged, while for a strong tube (run Ul, Bo = 8500 G), 67% of the tracers managed to 
emerge. This trend is due to the fact that the buoyancy of a flux tube scales with B 2 . 

The dependence of the percentage of emerged tracers on the twist parameter A is not so straightforward. 
We see that, at a given magnetic field strength of Bo = 8500 G, the percentage increases from 30% to 
68% in the range A = — 0.25. This trend is compatible with results from earlier idealized simulations of 
rising magnetic tubes. These studies found that an untwisted tube, after rising a distance of a fe w times its 
diam e ter, will split into counter-rotating vortex rolls which separate horizontally from each other (|Schusslerj . 
1 19771 ; lLongcope. Fisher, fc Arendtl Il996l ). The fragments of the individual tubes have smaller radii, less 
flux and are less buoyant. As a result, they are more easily disrupted by the external flow. Simulations 
of twiste d magnetic tubes show that the magnetic tension of the transverse field tends to keep the tube 



coherent dMoreno-Insertis fc Emonet . 1 19961; iFan. Zweibel. fc Lantzl . Il998l ; lEmonet fc Moreno-Insertisl . 1 19981 ; 



ICheung. Moreno-Insertis. fc Schiissleil . l2006t ). This is why the percentage of emerged tracers increases from 
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Figure 4: Same as Fig. but for run Ul (B n = 8500 G). In this run, the flux tube is sufficiently buoyant and 
the magnetic field sufficiently strong to rise against the drag of the downflow (see snapshot at t = min) . The 
expansion of the rising magnetic complex drives strong horizontal flows away from the emergence site. 
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Figure 5: Percentage of Lagrangian fluid elements initially embedded in the magnetic flux tube that reach 
base of the photosphere [z — 0) within the time interval < t < 20 min. 

A = to A = 0.25. For A > 0.25, however, the percentage drops again. This reversal in the trend is not in 
contradiction with the previous explanation. The complicating factor here is the relationship between the 
twist of a tube and its buoyancy. For the initial conditions we have chosen (i.e. the divergence of the total 
stress tensor at t = closely resembling that of the non-magnetic domain), the density deficit at the core of 
the tube scales as Ag/ (g) ~ 1/(1 + A 2 /2). Consequently, for a given initial field strength, the density deficit 
in the core of the tube actually decreases with increasing A. Hence, although the flux tube remains more 
coherent with increasing A, it is also somewhat less buoyant than the untwisted case. 

3.3 Trajectories of tracer fluid elements 

Figure [S] shows the vertical velocity and height of two tracers initially located at the axis of the magnetic 
flux tube. Solid and dashed lines correspond to results from run Ul and U6, respectively. The left panel 
corresponds to a tracer representing a segment of the initial flux tube that was embedded in a convective 
upflow. While the tracers in both runs accelerate and rise towards the surface, the tracer in run Ul (solid 
line) experiences larger acceleration and reaches higher rise speeds (v y « 2 — 3 km s _1 ) than the tracer in 
run U6 (v y ~ 1 km s" 1 ). This difference can be attributed to the higher buoyancy of the flux tube in run 
Ul. The right panel of Fig. |H]shows the corresponding results for a tracer in a tube segment initially located 
in a downdraft. Here, the tracer in run U6 is forced by the downdraft to descend by about 300 km before 
it gets caught up in a neighbouring upflow. The tracer in run Ul is also initially pushed downwards by the 
convective downdraft, but the flux tube in this case is sufficiently buoyant to counteract the downflow. 

Fig. [3 shows the time evolution of some further properties of the tracer initially embedded in an upflow 
(Fig. [6] left panel): plasma /3 (top panel), \B\/ g (middle), and Q ra d/f?, the radiative heating rate per unit 
mass. Until t ~ 8 min, the tracer fluid element rises almost adiabatically (Qrad ~ 0). During this time, 
as the fluid rises and expands, the value of /3 increases from an initial value of 1 to a maximum value of 
7. Thereafter, at about t = 8.5 min, the fluid element reaches the transparent photosphere and radiates 
intensely. The corresponding loss of internal energy corresponds to a decrease of the local gas pressure, 
leading to a compression of the fluid element and an abrupt decrease of /3 to a value of 3. Thereafter, the 
Lagrangian tracer continues to rise through the photosphere (eventually escaping through the top boundary 
of the domain). As it enters the more tenuous layers of the photosphere, it expands. Despite this, the plasma- 
/3 hovers for a few minutes at the level of /3 ~ 3 and then decreases even further to a value approaching unity. 
The time profile of the ratio \B\/g (middle panel of Fig. [7| provides a hint as to why this happens. We see 
that after the fluid has entered the photosphere, |-B|/<? begins to increase. By virtue of Walen's Equation 
(derived by combining the equation of continuity and the ideal induction equation, thus ignoring magnetic 
diffusion) , 



we see that |-B|/£> increases in the case of stretching of fluid along field lines and decreases in the case of 
de-stretching. 

We can compare the evolution of the tracer (which was initially embedded at the tube axis) with what 
would result from the simple model of a horizontal thin flux tube rising adiabatically in a stratified, static 




(9) 
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Figure 6: Motion of tracer fluid elements initially located at the axis of the magnetic flux tube. The left panel 
shows time profiles of rise velocity and height of a fluid element in a segment of the tube which was initially 
aligned with a convective upflow. The right panel shows the corresponding profiles for a fluid element in a tube 
segment in a convective downdraft. The solid and dashed lines indicate profiles for runs Ul and U6, respectively. 
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Figure 7: Temporal evolution of plasma /3, ratio \B\/g, and specific radiative heating, Q la d/ Q, for a tracer 
fluid element initially located at the axis of the tube segment embedded on a convective upflow (solid curves). 
The dashed curves indicates values corresponding to a thin flux tube model assuming adiabatic motion and 
no variation along the tube axis. As long as there is negligible radiative heating and draining of mass at the 
position of the tracer, there is a good match between the thin-tube model and the actual profiles. 
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atmosphere. The result is shown by the dashed curves in Fig. [7] The assumption that there are no variations 
along the direction of the tube axis means that \B\/g is constant (flat dashed line in the middle panel). 
Assuming adiabatic rise and instantaneous total pressure balance with the (prescribed) surroundings of the 
thin tube, we can determine the evo lution of the propert ies of the thin flux tube as functions of its initial 
state and its vertical displacement ()Cheung et all [2006). Between t = and t = 8 min, there is very 
good agreement between the actual and thin-tube model profiles, resulting from the quasi-adiabatic motion 
and negligible amount of mass drainage or accumulation along the field lines. Once the Lagrangian tracer 
reaches the surface and loses entropy, the profiles of (3 diverge. Whereas the actual value of (3 decreases 
sharply (due to cooling) and then more gradually (due to mass drainage) , the thin-tube model predicts that 
(3 continues to rise smoothly as the fluid continues to rise in the photosphere. The actual and model profiles 
of \B\/q also diverge. Since our thin-tube model does not take into account variations along the tube, \B\/ g 
remains constant in the model profile. In the actual profile, \B\/g increases after the fluid elements has 
emerged. This is ultimately a result of radiative cooling, which makes the emerged fluid denser than the 
fluid immediately below, driving the development of a downflow. This accelerating material in the developing 
downflow stretches the fluid along the field line and increases \B\/g. 



4 Observational signatures 



In the following sections, we compare diagnostics derived from the simulation results with observations of 
flux emergence. 

4.1 'Quiescent' flux emergence 

The tube in run U6 evolves passively with the convective flow. Consequently, the emergence of magnetic flux 
does not lead to a severe disturbance in the appearance of the granulation in white light intensity images. 
For this reason, we refer to these flux emergence events as quiescent. 

Figure[8]shows a sequence of continuum intensity and magnetic field maps representing the flux emergence 
event in run U6. The magnetic field is evaluated at T5000 =0.1 (optical depth of 0.1 for the visible continuum 
at 5000 A). We chose to sample from this optical depth because the maxima of the contribution functions 
of many photospheric lines commonly used to probe surface magnetic fields (e.g. Fe I 6301, Fe I 6302) 
are roughly located at this depth. As such, we can expect that the information carried by the line to be 
representative of the physical quantities in this layer of the atmosphere. In the left panels of Fig. [H] the 
overlaid colored contours indicate the vertical component of the magnetic field (B z ). 

The flux tube is initially aligned in the x-direction at y — 6 Mm. At t = 12.4 min, we see the first 
signs of flux emergence at the surface. The magnetic field emerges within the interior of granules, being 
predominantly horizontal near the granular centre. Towards the edge of the granule, the magnetic field 
becomes more vertical, indicating that the magnetic field threading a granule is forced to emerge in an 
arched configuration. One example is the bipolar region emerging in the granule centered at [x, y] = [4, 7] 
Mm (see snapshot at t — 12.4 min). At this particular instant, the flux contained in each polarity of this small 
bipole is 8 x 10 17 Mx. The diverging horizontal flow of the granule flow has the effect of expelling the emerged 
flux in to the network of intergranular dow nflow lanes within a typical granulation timescale on the order of 



5 min (jBrav et al ], Il984l : iTitle et al.l,ll989l). This finding is in accordance with earlier simulations o f magneto- 



convection in an idealized setup (|Weissl . 1 19661 : IProctor fe Weissl . [l982 ; Hurlb urt fe Toomrd . [ 1988) as well as 



more recent, realistic simulations of magneto-convection i n the near-surface layers of the solar conve ction 
zone ijSteiner et allll998l : lB~ercikl . l2002l : [Bercik et alll2003l : FvHgler et al.l . [2005l : IStein fe NordlundLl2006h . At 
t — 17.7 min (about 5 minutes after the first snapshot), the distribution of the vertical field B z already 
shows that the regions of strongest field are in the intergranular lanes. However, the surface field strength 
at T5000 =0.1 rarely exceeds the equipartition value of 450 G. In comparison, the granule interiors have 
significantly weaker fields. The next snapshot at t = 22.1 min illustrates the further transport and dispersion 

of the emerged flux in the intergranular network. 

How do our simulation results compare with observations? iDe Pontieul (|2002l ) has obtained high cadence 
and high resolution observations of an emerging flux region (EFR). Within a time series of 3.5 hr duration, 
he identified seven individual small-scale emergence events. His two key findings relevant for this discussion 
are: (a) magnetic concentrations emerge in the interior of granules; and (b) within 10 — 15 min from their 
initial appearance, the flux concentrations are quickly dispersed by the granular flow. He estimates the 
average flux density of the emerging magnetic concentrations to be about 200 ± 30 Mx cm -1 . If this value 
were to correspond to the actual magnetic field strength, the emerging fiel ds would ha ve sub-equipartition 
strength (B cq ~ 450 G), a finding consistent with his results (a) and (b). IDe Pontieul also gives estimates 
for the average flux of the emerging concentrations in the range (9 ± 4) x 10 17 Mx. This is the same order 
of magnitude as what we find for individual polarities of emerging flux concentrations in run U6. 
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Figure 8: Emergence of magnetic field into the photosphere in run U6. Left column: time sequence of the 
vertically emerging continuum intensity at 5000 A, Isooo- Right column: the corresponding time sequence 
of magnetic field inclination evaluated at T5000 =0.1. A field inclination (also called zenith angle) of 0° 
and 180° respectively correspond to vertical magnetic field directed out of and into the x — y plane. In the 
left panels, the color contours indicate the vertical component of the magnetic field at T5000 =0.1 with the 
levels ±[50,100,200,400] G. Green (red) contours indicate magnetic upward (downward) directed field. An 
accompanying mpeg animation is available at http://www.mps.mpg.de/homes/cheung/U6emergence.mpg, 
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Figure 9: Left column: Joint probability distribution functions (JPDFs) of the zenith angle of the magnetic 
field vector and the absolute field strength, \B\, both at the level T5000 = 0.1. Right column: JPDF of zenith 
angle and vertical fluid velocity. The colour coding indicates the logarithm of the probability. The upper and 
lower rows, respectively, show the JPDFs at t = 12.4 min and t = 61.1 min for run U6. 
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|Pe Pontieul suggested that the small-sc ale emergence eve nts he has detected may be related to the so- 
called horizontal internetwork fields (HIFs. lLites et aD . 1 19961 s ! . which are predominantly horizontal magnetic 
structures with angular and time scales of 1 — 2" and 5 min, respectively. The existence of HIFs are inferred 
from the distinctive signals in the Stokes Q and U profiles of relatively quiet regions of the solar surface. 
The strength of these signals indicate that HIFs are generally weak with \B\ < 600 G. HIFs are usually 
detec ted in regio ns w here the spe ctral lines are slightly blueshifted, indicating an association with upflows. 
Both iLites et al.l and |Pe Pontieul interpret these as signatures of the crests of magnetic loops emerging in 
granules. 

We have analyzed simulation run U6 to check whether this interpretation is consistent with our results. 
Fig-E]shows joint probability distribution functions (JPDFs) of different surface quantities in the simulation 
run for times t — 12.4 min (upper row) and t — 61.1 min (lower row). The left column shows JPDFs between 
the zenith angle 7 of the magnetic field and the absolute field strength, \B\, both evaluated at T5000 = 0.1. 
The zenith angle is defined as the angle between B and the vertical direction (z-axis). As such, magnetic 
field vectors with 7 = 0° and 7 = 180° correspond to purely vertical magnetic field that point out of and into 
the horizontal plane, respectively. Purely horizontal fields have 7 = 90° . The right column shows JPDFs of 
the zenith angle of B and the vertical fluid velocity (for surface regions with \B\ > 0). Flux begins to emerge 
approximately at t — 12.4 min. The JPDFs at this time indicate the existence of predominantly horizontal 
field with strengths of up to \B\ ~ 400 G. These horizontal fields mainly have rise velocities of 1 — 2 km s -1 , 
which is typic al of granular up flows. This finding is compatible with the scenario proposed bv lLites et alj 
l|l996h and bv lDe Pontieul l|2002ft for the explanation of HIFs. 

The snapshot taken at t = 61.1 min (lower panels in Fig. [9] shows a distinctly different shape for the 
JPDF between the zenith angle and \B\. The pair of horn-like features indicate strong vertical fields of up to 
1 kG. We will return to this point in the following section. Furthermore, there is a continuous distribution 
of nearly horizontal fields with strengths of up to 300 G. The accompanying JPDF between the zenith angle 
and v z indicates that the majority of horizontal fields reside within upflow regions (granules). In fact, the 
JPDF between the zenith angle and \B\ at t = 61.1 min resembles those found in simulations of magneto- 
convection in the near-surface layers with a net vertical magnetic flux threading the simulation domain 
l|Voglerl . 12003? ) . This indicates that the abundance of weak horizontal fields in granule interiors is a general 
feature of magneto-convection and not necessarily always associated with emerging flux. 

4.1.1 Surface evolution of emerged field: cancellation, coalescence and secondary 
emergence 

In run U6, the morphology of the emerged field bears little resemblance to the initial horizontal flux tube 
structure. The arrangement of the surface field after the initial flux emergence (see snapshots at t = 17.7 
an d 22.1 min in Fig. El) resembles th e 'salt and pepper' pattern of quiet Sun magnetic fields observed, e.g., 
by iDomfnguez Cerdena et al.l (|2003h ■ displaying a mixture of positive and negative small-scale flux in the 
intergranular network. These authors describe instances when flux concentrations of like polarities meet 
and appear to coalesce (to within the resolution limit). They also show a flux cancellation event due to 
the encounter of opposite polarities. In our simulation run U6, both these types of events occur as well. In 
the following, we discuss in detail an instance of a coalescence event between vertical flux concentrations of 
like-polarity. 

Figure [TU] shows a time sequence of intensity images in the vicinity of a granule, overplotted with hor- 
izontal velocity vectors and contours of the vertical magnetic field strength, both at T5000 = 0.1. At time 
t — 47.9 min, we find a number of magnetic flux concentrations of both polarities. The most prominent 
feature is centered at [x,y] = [2.5,9] Mm. This feature resides at a downfiow vertex of intergranular lanes 
and has a flux of 1.5 x 10 18 Mx. The core of the feature has a maximum vertical field strength of B z = 1, 200 
G (at T5000 = 0.1). A few minutes later, at t = 52.1 min, we find a small bipole emerging through a granule 
that is located to the right of the pre-existing flux concentration. The bipole is oriented such that the pole 
closest to the pre-existing flux concentration has the same polarity (positive). The newly emerged flux is 
subsequently expelled to the adjacent downfiow lane, where it coalesces with the pre-existing feature. The 
resulting larger flux concentration has a net flux of 2.6 x 10 18 Mx, comparable to the longitudinal flux of 
the horizontal flux tube at the beginning of the simulation ($0 = 3.1 x 10 18 Mx). 

This particular example is interesting for a number of reasons. Firstly, the new bipole emerges at a 
relatively late stage of the simulation, about 40 min after the initial appearance of flux at the surface (see 
Fig. [8} . Its emergence location is also unusual. In Fig. [8] we see that almost all the flux emerges between 
y = 4 and y = 8 Mm. The new bipole, however, emerges significantly further away, at y — 9.2 Mm. These 
two points are consistent with the fact that, when a section of a flux tube is passively carried by an upwelling, 
not all the flux contained in the tube emerges in one single event. The material near the edge of an upwelling 
may overturn before it reaches the surface. Consequently, the magnetic field threading this material fails to 
emerge on first attempt. This magnetic field may then continue to travel downwards, or may get caught up 
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Figure 10: Coalescence of newly emerged magnetic concentrations with pre-existing surface field in simulation 
run U6. The grey shading in each panel represents the normalized continuum intensity at 5000 A. The vector 
overlays indicate the horizontal velocity field at optical depth T5000 = 0.1 and the contour overlays indicate B z 
at the same optical depth at the levels ±[100, 200, 400, 800, 1600] G. Green (red) contours correspond to positive 
(negative) polarities. At t — 52.1 min, a bipole centered at [x, y] — [3.5, 9.2] Mm is emerging through a granule. 
Flux expulsion leads to coalescence of the pre-existing positive flux with the positive part of the newly emerged 
flux. 



in another upwelling and travel upwards again. The emergence of the bipole in Fig. [T0]is one such example. 
The magnetic field associated with this bipole has overturned several times before it eventually reaches the 
T5000 =0.1 surface. Considering that it had been caught up in different upflows and downflows, it is not 
surprising that its emergence location is about 1 Mm further afield. 

The coalescence event shown in Fig. \W\ is interesting for another reason. The pre-existing positive 
flux concentration residing in a downflow vertex is not as dark as the average intergranular network. After 
coalescence, the resulting concentration becomes even brighter. To explain this, we make reference to Fig. 1111 
Let us first consider the top row of this figure, which shows joint probability distribution functions (JPDFs) 
between the value of \B Z \ at a horizontal plane of constant geometrical height (zo.i) and its value evaluated 
at optical depth T5000 = 0.1. The height (20.1) is defined as the average geometric height of the T5000 = 0.1 
surface in the absence of magnetic fields. On average, the T5000 = 0.1 surface lies about 160 km above the 
level T5000 = 1.0. We first focus on the distribution of field strengths at the horizontal plane z = (20.1) 
(i.e. ignore the y-axis). A few minutes after the initial appearance of flux at the surface (t = 17.7 min), 
the vertical field strengths do not exceed 500 G. Some time later, at t = 34.6 min, a small fraction of the 
vertical flux concentrations have reached field strengths of up to 700 G, and at t = 61.1 min we find fields 
approaching \B Z \ = 1 kG at z = (20.1). 

The appearance of stronger magnetic fields several g ranulation timescales after the initial emergence of 
flux suggests that c o nvect i ve intensification operates (e.g.. lNordfund fc Steinlll990l ; [Grossmann-Doerth et all 
1 19981 ; IVogler et all 120051 : IStein fc Nordlundt |2006j h The magnetic field maps from run U6 show that the 
maximum value of \B Z \ within a vertical magnetic bundle is an increasing function of its total magnetic 
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Figure 11: Joint probability distribution functions (JPDFs) for various quantities in run U6. In all panels, the 
horizontal axis gives the value of \B Z \ at a plane of constant geometrical height, z — (zo.i)- The quantities given 
by the vertical axis are: \B Z \ at the local level T5000 = 0.1 (first row of panels), the local vertical displacement 
of the T5000 = 0.1 surface from the mean geometrical height of the surface corresponding to this optical depth 
( z o. 1 — ( z o.i), second row), the temperature at T5000 = 1-0 (third row), and the normalized emergent continuum 
intensity at 5000 A(fourth row). The three columns show the JPDFs at three different times. The color coding 
indicates the logarithm of the probability. 
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Figure 12: Scatter plot of the maximum value of \B Z \ in individual magnetic flux concentrations versus the 
amount of vertical flux contained in the flux concentration. The values of \B Z \ were evaluated at a plane of 
constant geometrical height, z — (zo.i). The crosses, diamonds and squares indicate values corresponding to 
the snapshots taken at t = 17.7 min, t = 34.6 min, and t = 61.1 min, respectively. 



flux, <F In order to identify individual flux bundles, we choose a threshold value, B t h. For each discrete 
region in the magnetic field map with \B Z \ > B t h, we measured the maximum value of \B Z \ in the region as 
well as the amount of flux contained within the region. A scatter plot of these two quantities is shown in 
Fig-EU The crosses, diamonds and squares correspond to the snapshots taken at t — 17.7 min, t = 34.6 min, 
and t = 61.1 min, respectively, for B t h — 50 G. The scatter plot clearly shows that larger flux bundles can 
support stronger internal field s trengths. This trend is n ot sensitive to the threshold value used. Our finding 
is compatible with the work of IVenkatakrishnanl (| 1986T ) - who suggested that the radiative heating of flux 
tubes by their surround i ngs in hibits the intensification of very small flux tubes - and with the observational 
study bv ISolanki et al.1 |l996:). We should point out, however, that the various types of diffusion present 
in the simulation affect the results for the smallest flux concentrations. For instance, a flux bundle with a 
mean field strength of 100 G and a flux of 10 16 Mx has a radius of about 50 km, equal to the horizontal grid 
spacing used in our simulations, so that it is essentially unresolved. Values in the scatter plot in the range 
$ < 10 17 Mx are therefore probably influenced by diffusive effects. 

For the snapshot taken at t — 61.1 min, the maximum field strength at z = (zo.i) has reached about 
1 kG. The corresponding magnetic pressure is almost equal to the horizontally averaged gas pressure at that 
geometrical level. This indicates that the strongest flux concentrations are largely evacuated, leading to a 
downward shift of the optical depth scale within the flux concentration ('Wilson depression'). The relation 
between the displacement, Zo.i — (zo.l), and the field strength is shown in the second row of JPDFs in Fig. 1111 
the stronger the field, the larger is the downward (negative) displacement of the Tgooo =0.1 level, which can 
reach values up to 120 km. This has the consequence that \B Z \ at the local level T5000 = 0.1 has values of 
up to 1.5 kG, since the field strength increases with geometrical depth (see the top right panel in Fig. Hip . 

The partial evacuation of vertical flux concentrations leads to another observational signature. From 
the Eddington-Barbier relation, we know that the emergent intensity is approximately given by the source 
function at r = 1. The panels in the third row of Fig. 1111 show JPDFs between \B Z \ at z = (zo.i) and the 
temperature at T5000 = 1-0. The temperature at optical depth unity is higher in strong field regions, with the 
result that the brightness of vertical flux tubes increases with \B Z \ (see fourth row of Fig. Hip . The higher 
temperature at T5000 = 1.0 i s due to the l ateral radiativ e heating of the interior of the flux concentration by 
the s urrounding 'hot walls ' l|Spruitl[l976l ; lDeinzer et all 1 1984 Ig rossmann- Doerth et allll994lVogi er et al.1 . 
l2005h . 



4.2 Emergence of strong magnetic field 

In this section, we discuss the observational signatures of the emerging flux tube with Bo ~ 8500 G (run Ul), 
which is strong enough to resist being passively deformed and advected by the convective flow. Fig. |13(a)| 
gives a time sequence of continuum intensity maps throughout the emergence event (grey scale). The 
vector overlays indicate the horizontal components of the velocity field at T5000 = 0.1. Fig. 13(b) shows the 
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Figure 13: Flux emergence event in run Ul. Left: an arrow with a length of one grid-spacing on the vector plot 
corresponds to a horizontal flow speed of v = 3 km s -1 . Right: an arrow with a length of one grid-spacing on 
the vector plot corresponds to a horizontal field strength of \B\ — 400 G. Both the velocity and magnetic fields 
shown are evaluated at T5000 = 0.1. The granulation pattern is temporarily disturbed by the emerging magnetic 
field and a transient darkening along the length of the domain marks the site of emergence. An accompanying 



mpeg animation is available at |http://www.mpsTm pg.de/homes/cheung/Ulemergence.mpg 



corresponding maps of the vector magnetic field. 

We have seen that the emergence morphology of the field in run U6 is dictated by the convective motion 
of the near-surface layers. In run Ul (where the magnetic field is several times stronger), the convection 
still has an appreciable influence on the emergence morphology. Comparison of the magnetic field and the 
intensity maps at t = 8.6 min clearly shows that magnetic flux emerges preferentially within the interior of 
granules. On the other hand, we find that the emerging magnetic structure also affects the appearance of 
the local granulation pattern. At t = 8.6 min (top panel of Fig. 13(a) I, there is an alignment of granules 
along the direction of the tube axis. These granules are relatively large, each covering an area of up to 9 
Mm 2 . The disturbance of the granulation by the emerging magnetic field is most prominent at t — 11.2 
min, showing a darkening near the emergence site extending over the full horizontal size of the simulation 
domain. This darkening is not coherent along its entire length: there are bright patches on either side of 
the darkening but their shapes are distinctly different from those of normal granules. The expansion of the 
emerging tube drives horizontal outflows away from the emergence site, which are visible in the vector plots 



in the middle panel of Fig. 13(a) The magnetic field maps for the same instant shows that much of the 
emerged flux already resides in intergranular lanes. 

At t = 16.9 min, the darkening has largely vanished and a granulation pattern more similar to the normal 
pattern has re-established itself. However, these granules tend to be elongated in the y-direction because of 
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Figure 14: Joint probability distribution functions (JPDFs) of the zenith angle of the magnetic field vector 
versus the absolute field strength \B\ in run Ul. Both quantities have been evaluated at the T5000 = 0.1 level. A 
zenith angle of 90° corresponds to a purely horizontal field. The four panels show the JPDFs at four different 
times. At t — 8.6 min, the flux tube is emerging at the surface and the strongest fields are predominantly 
horizontal. Kilogauss fields develop only later and are almost exclusively vertical at 27.6 min. 

the horizontal outflows driven by the expanding tube. In the corresponding magnetic field map, we find that 
the emerged magnetic flux outlines the boundaries of these new granules. The y-component of the magnetic 
field is, on average, positive. This is just the opposite to what we see in the magnetogram at t = 8.6 min. 
This reversal of the sign of B y within this time span is expected. In the earlier snapshot, the T5000 =0.1 
surface intersects with the upper half of the flux tube. The twist of the field lines within the tube is such 
that By < in this part of the tube. In the latter snapshot, the lower half of the tube has emerged at the 
surface, and so the synthetic magnetogram shows field with predominantly B v > 0. 

4.2.1 The relation between field strength and zenith angle 

1 1 j k 

ILites. Skumanich. fc Martinez Pilletl ([1998) have presented observational results for several emerging flux 
regions. They obtained vector magnetogra ms by applying in version methods to measurements of the full 
Stokes vector. Based on their observations. ILites et al.l (|l998h describe the following scenario for flux emer- 
gence: magnetic flux emerges at the surface in the form of horizontal field with field strengths of 200 — 600 G. 
After emergence, the field quickly migrates away from the emergence site. In the process, the field becomes 
ve rtical. Only then d o the fields obtain field strengths exceeding one kilogauss. The observational study 
bv lKubo et all (|2003h also supports these findings. 

Figure [14] shows JPDFs of the zenith angle of the magnetic field vector versus the absolute field strength 
from simulation run Ul. At t = 8.6 min, the flux tube is just emerging at the surface. The JPDF at this time 
reflects the fact that the flux tube is still coherent. The horizontal fields are stronger than the vertical fields. 
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Figure 15: Dark lane (y ~ 6 Mm) in run Ul ( £ = 11.2 min). The grey shading represents the normalized emer- 
gent continuum intensity. The blue contours indicate upflows with v z = 0.5 km s _1 and 1 km s _1 , respectively, 
while the red contours correspond to downflows of the same magnitude (v z < 0). The dark lane is roughly 
coincident with upflowing material. 



Although a tiny fraction of the horizontal fields have kilogauss field strengths, most are confined within the 
range 400 < |S| < 1000 G. Just five minutes later, at t — 11.2 min, the distribution in the JPDF looks 
very different. The appearance of two 'horns' in the JPDF near 7 = 0° and 7 = 180° indicates the effect 
of the granular flow on the emerged field. As the emerged field is advected to the boundaries of granules, 
it also becomes predominantly vertical. As the simulation progresses in time, the two horns in the JPDF 
become more distinct. At t = 27.6 min, nearly all the field with \B\ > 600 G is vertical. There is actually 
an abundance of horizontal fields, but these are confined to strengths of less than 6 00 G. The shape o f the 
JP DF at this t i me is similar to the scatter plots of zenith angle vs. |S| obtained bv iLites et al l (|l998h and 
bv lKubo et all (|2003l h 



4.2.2 The dark lane 



The physical nature of the extended dark lane (see central panel of Fig. 13(a) I is distinctly different from 
that of normal intergranular lanes, which are associated with cold downflowing material. Instead, the dark 
lane that results from the emergence of the flux tube in run Ul is associated with upflowing material for 
a substantial fraction of its lifetime (several minutes) . Fig. [TS] shows the normalized emergent continuum 
intensity at t = 11.2 min (same as the central panel in Fig. I10[l . Overplotted on the grey scale image 
are contours of the vertical component of the velocity at T5000 = 0.1 for magnetic regions (\B\ > 400 G). 
The blue contours indicate upflows while the red contours correspond to downflows. There is a substantial 
overlap of the blue contours with regions on the surface comprising the dark lane. The correlation between 
upflows and low emergent intensity in the dark lane is also illustrated in Fig. 1161 which shows profiles of 
emergent continuum intensity and vertical velocity at T5000 = 0.1 across the dark lane (in the y-direction). 
The plotted values are averages taken along the :r-direction at time t — 11.2 min. Away from the emergence 
site, we have the usual intensity-velocity correlation of granules and intergranular lanes. However, within 
the emergence site (around y ~ 6 Mm), the darkening (/sooo/(^500o) ~ 0.93) is co-spatial with upflowing 
material (v z — 0.5 — 1 km s" 1 ). 

The lifetime of the extended dark lane is about 10 minutes, which is comparable to the timescale of 
granulation. During the first few minutes, the regions of lower brightness are associated with upflowing 
material with rise speeds of up to 1 km s — 1 . This is substantially lower than the typical rise speed (v ~ 
2 — 4 km s _1 ) of magnetic material when it is just emerging at the visible surface. In the last few minutes of 
the dark lane's lifetime, the dark lane splits up into spatially separated, dark elongations with length of about 



2 Mm (see intensity image at t = 16.9 min in Fig. 13(a) I. The vertical velocities of these dark elongations are 



negative, i.e., they have developed into downflows. The presence of the dark lane depends on the initial field 
strength as well as the initial twist carried by the tube. When the initial flux tube is sufficiently weak that 
its evolution is dominated by the drag force (as is the case with run U6), a dark lane does not appear: the 
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Figure 16: Average profiles of the normalized continuum emergent intensity (solid line, left axis) and vertical 
velocity at r 50 oo = 0.1 (dashed line, right axis) across the dark lane at t = 11.2 min. The curves represent 
averages over the x-direction. 



magnetic flux is passively carried upward by the convective flow, which shows the normal intensity-velocity 
correlation. In order for the dark lane to appear, the emerging field must be strong enough that its evolution 
is not completely controlled by the convective flow. However, this is only a necessary criterion. The amount 
of initial twist is also important. We have carried out simulation runs with the same initial conditions for 
the flux tube as in run Ul, varying only the amount of the initial twist. We found that for a flux tube with 
twist parameter A = 0, there is no transient dark lane marking the emergence site of the tube. For a flux 
tube with a moderate amount of twist, A = 0.25, we find the transient appearance of some dark patches 
within the granules where the emergence occurs. These patches, however, do not exhibit a clear alignment 
or coherence as is the case in run U6 (where A = 0.5). The following explanation accounts for the properties 
we have presented above: When a flux tube (with a flux of, say, $o = 10 19 Mx) reaches the surface, it 
has a rise speed of 2 — 4 km s _1 . As demonstrated by previous numeric al MHD simulations of the rise of 
twisted flux tubes through an init i ally static atmosphere (see, for example .iMoreno-Insertis fc Emonetl 1 19961 : 
lEmonet fc Moreno- Insertisl. Il998l ; ICheung. Moreno-Insertis. fc Schusslej . l2006l ). the amount of twist in the 
tube determines whether it can remain a coherent structure (large A) or break up into vortex fragments 
(small A). As the flux tube emerges and cools by radiation, it loses buoyancy. Nevertheless, the emerged 
magnetic material may still overshoot provided that the flux tube is sufficiently twisted, so that its core 
remains largely intact and maintains its magnetic buoyancy. As the emerged magnetic material of the tube 
overshoots, its upwards motion is braked by the stable stratification. This is why the upward velocities in 
Fig. are limited to 0.5 — 1 km s _1 within the dark lane. The strong cooling of the magnetic material 
within the core of the tube at optical depth unity accounts for the lane's darkness. The arrangement of 
the cool and dense material in an extended lane configuration is not stable to perturbations introduced by 
granular dynamics. Within a few minutes, the extended dark lane breaks up into shorter elongations with 
lengths comparable to those of intergranular lanes. The buoyancy loss of the material by radiation and 
stratification either leads to the birth of new downflows, or feeds pre-existing downflows in the intergranular 
network. The influence of the granular dynamics on its evolution explains why the dark lane has a lifetime 
comparable to the granulation timescale. 

Transient appe arance of dark alignments in emerging flux regions have been reported by a number of ob- 
servational studies dBrav. Loughhead fc Durrantl ll984l ; lBrants fe Steenbeekl ll985l : IZwaar] , Il985l ; IStrous fe Zw aan , 
Il999l h For example. IStrous fe Zwaanl l|l999^ have carried out a statistical analysis of an emerging flux region 
and found that transient darkenings in the continuum as well as in the cores of photospheric spectral lines, 
lasting about 10 minutes, are robust signatures of emergence events. They note that the typical extension 
of the darkenings is about 2 Mm and that they are aligned ro ughly parallel to t he dir ection connecting the 
two developing polarities of the corresponding active region. IStrous fc Zwaanl ljl999h also report that the 
darkenings are associated with upflows of about 0.5 km s -1 . T hese properties are reproduced by our simu- 
lation results from run U6. In addition. IStrous"fe Zwaanl (|l999T ) also report that the darkenings, on average, 
rotate counterclockwise at about 0.5° min -1 . In run U6, we do find that some of the dark elongations rotate 
during the breakup of the dark lane. However, we do not find any systematic direction in which the rotation 
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occurs. 



5 Discussion 

The magnetic flux emergence simulations presented in this paper are 'realistic' in two senses. Firstly, 
the simulations take into account the effects of non-local radiative energy exchange, partial ionization and 
magneto-convection, all of which are important for a proper treatment of the problem. Secondly, the 
simulations yield observational signatures of magnetic flux emergence that are in qualitative and quantitative 
agreement with observations of emerging flux. 

Our simulation results are generally consistent with the findings of iFan. Abbett. fc Fisherl (|2003h . ex- 
tending their validity to the near-surface layers, where the anelastic approximation used by these authors is 
no longer tenable. However, Mach numbers reaching M ~ O(0. 1) — 0(1) and the steep decrease of pressure 
and density with height limit the range of possible ratios of initial and equipartition field strengths, Bo/B cq . 
As a consequence, a considerable influence of the convective flows on the dynamics of the rising flux tube 
and on the flux emergence process is unavoidable. 

Convection influences the evolution of emerging magnetic fields both before and after they appear at the 
photosphere. In the subphotospheric layers, the convective flow distorts a rising flux tube and imposes a 
systematic undulation along the tube (with a typical wavelength ~ 1 — 2 Mm). Whereas upflows aid the 
rise of certain segments of the tube, downflows supp ress the rise of othe r segments. Our work shows that 
the undulatory nature of an emerging magnetic tube (IPariat et all [2004) is a consequence of its interaction 
with granular convection. 

The amount of distortion and deformation suffered by a magnetic flux tube interacting with convective 
flows depends on its initial properties, most notably its buoyancy and twist. Magnetic flux tubes with less 
longitudinal flux (say, 10 18 Mx) are not sufficiently buoyant to rise coherently against the convective flows. 
The emergence events associated with these smaller and weaker flux tubes are inconspicuous, in the sense 
that the granulation pattern in the quiet Sun is not disturbed (see Section [TTTJ . Their small spatial scales 
(~ 1 Mm) and short tempo ral scales (5 min) make their observations difficult. Although detections of such 
events have been reported l|De PontieuLliooa i. it is unknown how much flux emerges at the surface in this 
form. High-cadence, high-resolution observations from the Solar Optical Telescope onboard Hinode (Solar-B, 
llchimoto et all 120041 ; iTarbelll . 120061 ') will allow us to study such events in more detail then ever before. 

In Section 14.1.11 we provided an example of the secondary emergence of a bipole. This case (with a flux 
of 10 18 Mx in each polarity) is peculiar for two reasons. Firstly, the bipole emerged rather distant from 
the main emergence site. Secondly, it emerged several granulation time scales after the initial appearance 
of flux at the surface. This secondary emergence event is the result of the recirculation and overturning of 
material in the near-surface layers of the convection zone. It suggests to us that part of the small-scale flux 
emergence events on the solar surface could result from recirculation of material in the convection zone. 

The expulsion of emerged magnetic flux from the granule interiors to the intergranular network on the 
order of the granulation time scale (5 — 10 min) is a common feature of all our flux emergence simulations. 
Having said this, not all emergence events appear the same. For example, the observational signatures of 
the emergence of tubes with 10 19 Mx flux are clearly different from those of tubes with 10 18 Mx. Although 
they still suffer distortion and deformation by the convective flow, such bigger flux tubes are able to rise 
against downflows and eventually emerge at the surface. During emergence, the strong horizontal expansion 
of rising fluid can disturbs the granulation pattern. If the flux tube is sufficiently twisted (say, A = 0.5), the 
emergence event is accompanied by the transient appearance of an extended darkening in the direction of 
the tube axis. Such features have been reported in observations of emerging flux regions. Analysis of our 
simulations results reveals that they are associated with emerged, cooled magnetic material overshooting 
into the photosphere. 

Acknowledgements 

The authors thank the anonymous referee for helping to substantially improve the presentation of this paper. 
MCMC wishes to thank all the friends and colleagues at the IAC for their hospitality during his extended 
stays in Tenerife as a visiting graduate student. In particular, MCMC is deeply grateful to have known the 
late Juan Luis Medina Trujillo, who will always be remembered for his sincerity, his love of nature and his 
zest for life. 

References 

Abbett, W. P. & Fisher, G. H. 2003, ApJ, 582, 475 
Abbett, W. P., Fisher, G. H., & Fan, Y. 2000, ApJ, 540, 548 



23 



Archontis, V., Moreno-Insertis, P., Galsgaard, K., Hood, A., & O'Shea, E. 2004, ApJ, 441, 886 
Archontis, V., Moreno-Insertis, F., Galsgaard, K., & Hood, A. W. 2005, ApJ, 635, 1299 
Bercik, D. J. 2002, PhD thesis, AA(Michigan State University) 

Bercik, D. J., Nordlund, A., & Stein, R. F. 2003, in ESA SP-517: GONG+ 2002. Local and Global Helio- 
seismology: the Present and Future, ed. H. Sawaya-Lacoste (ESA Publications Division), 201-206 

Born, R. 1974, Sol. Phys., 38, 127 

Brants, J. J. & Steenbeek, J. C. M. 1985, Sol. Phys., 96, 229 

Bray, R. J., Loughhead, R. E., & Durrant, C. J. 1984, The solar granulation (2nd edition) (Cambridge and 
New York, Cambridge University Press, 1984, 270 p.) 

Caligari, P., Moreno-Insertis, F., & Schiissler, M. 1995, ApJ, 441, 886 

Caligari, P., Schiissler, M., & Moreno-Insertis, F. 1998, ApJ, 502, 481 

Cheung, M. 2006, PhD thesis, University of Gottingen, Germany, 

http://webdoc.sub.gwdg.de/diss/2006/cheung 

Cheung, M. C. M., Moreno-Insertis, F., & Schiissler, M. 2006, A&A, 451, 303 

Cheung, M. C. M., Schiissler, M., & Moreno-Insertis, F. 2007, A&A, 461, 1163 

De Pontieu, B. 2002, ApJ, 569, 474 

Deinzer, W., Hensler, G., Schiissler, M., & Weisshaar, E. 1984, A&A, 139, 426 
Dommguez Cerdena, I., Sanchez Almeida, J., & Kneer, F. 2003, A&A, 407, 741 
Dorch, S. B. F., Gudiksen, B. V., Abbett, W. P., & Nordlund, A. 2001, A&A, 380, 734 
D'Silva, S. & Choudhuri, A. R. 1993, A&A, 272, 621 
Emonet, T. & Moreno-Insertis, F. 1998, ApJ, 492, 804 
Fan, Y. 2001, ApJ, 554, Llll 

Fan, Y., Abbett, W. P., & Fisher, G. H. 2003, ApJ, 582, 1206 
Fan, Y., Fisher, G. H., & Deluca, E. E. 1993, ApJ, 405, 390 
Fan, Y., Fisher, G. H., & McClymont, A. N. 1994, ApJ, 436, 907 
Fan, Y. & Gibson, S. E. 2004, ApJ, 609, 1123 
Fan, Y., Zweibel, E. G., & Lantz, S. R. 1998, ApJ, 493, 480 
Forbes, T. G. & Priest, E. R. 1984, Sol. Phys., 94, 315 

Galsgaard, K., Archontis, V., Moreno-Insertis, F., & Hood, A. 2007, ApJ, submitted 
Galsgaard, K., Moreno-Insertis, F., Archontis, V., & Hood, A. 2005, ApJ, 618, L153 
Grossmann-Doerth, U., Knolker, M., Schiissler, M., & Solanki, S. K. 1994, A&A, 285, 648 
Grossmann-Doerth, U., Schiissler, M., & Steiner, O. 1998, A&A, 337, 928 
Hagenaar, H. J. 2001, ApJ, 555, 448 

Hagenaar, H. J., Schrijver, C. J., & Title, A. M. 2003, ApJ, 584, 1107 
Hurlburt, N. E. & Toomre, J. 1988, ApJ, 327, 920 

Ichimoto, K., Tsuneta, S., Suematsu, Y., et al. 2004, in Optical, Infrared, and Millimeter Space Telescopes. 
Edited by Mather, John C. Proceedings of the SPIE, Volume 5487, pp. 1142-1151 (2004)., ed. J. C. 
Mather, 1142-1151 



24 



Isobe, H., Miyagoshi, T., Shibata, K., & Yokoyama, T. 2005, Nature, 434, 478 
Kubo, M., Shimizu, T., & Lites, B. W. 2003, ApJ, 595, 465 
Leake, J. E. & Arber, T. D. 2006, A&A, 450, 805 

Lites, B. W., Leka, K. D., Skumanich, A., Martinez Pillet, V., & Shimizu, T. 1996, ApJ, 460, 1019 

Lites, B. W., Skumanich, A., & Martinez Pillet, V. 1998, A&A, 333, 1053 

Longcope, D., Fisher, G., & Arendt, S. 1996, ApJ, 464, 999 

Magara, T. 2001, ApJ, 549, 608 

Magara, T. 2006, ApJ, 653, 1499 

Magara, T. & Longcope, D. W. 2003, ApJ, 586, 630 

Moreno-Insertis, F. 1983, A&A, 122, 241 

Moreno-Insertis, F. & Emonet, T. 1996, ApJ, 472, L53 

Murray, M. J., Hood, A. W., Moreno-Insertis, F., Galsgaard, K., & Archontis, V. 2006, A&A, 460, 909 

Nordlund, A. & Stein, R. F. 1990, in Solar Photosphere: Structure, Convection and Magnetic Fields, IAU 
Symp. 138, ed. J. O. Stenflo (Dordrecht: Kluwer Academic Publishers), 191 

Pariat, E., Aulanier, G., Schmieder, B., et al. 2004, ApJ, 614, 1099 

Parker, E. N. 1955, ApJ, 121, 491 

Parker, E. N. 1979, ApJ, 230, 914 

Proctor, M. R. E. & Weiss, N. O. 1982, Reports of Progress in Physics, 45, 1317 
Roberts, B. & Webb, A. R. 1978, Sol. Phys., 56, 5 
Schiissler, M. 1977, A&A, 56, 439 
Schiissler, M. 1979, A&A, 71, 79 

Schiissler, M. 1984, in The Hydromagnetics of the Sun (European Space Agency, ESA SP-220), 67-76 

Schiissler, M. 1987, in ASSL Vol. 137: The Internal Solar Angular Velocity, ed. B. R. Durney & S. Sofia, 
303-320 

Shibata, K., Tajima, T., Steinolfson, R. S., & Matsumoto, R. 1989, ApJ, 345, 584 

Solanki, S. K., Zufferey, D., Lin, H., Rueedi, I., & Kuhn, J. R. 1996, A&A, 310, L33 

Spruit, H. C. 1976, Sol. Phys., 50, 269 

Spruit, H. C. 1981, A&A, 98, 155 

Stein, R. F. & Nordlund, A. 1989, ApJ, 342, L95 

Stein, R. F. & Nordlund, A. 1998, ApJ, 499, 914 

Stein, R. F. & Nordlund, A. 2006, ApJ, 642, 1246 

Steiner, O., Grossmann-Doerth, U., Knolker, M., & Schiissler, M. 1998, ApJ, 495, 468 

Strous, L. H. & Zwaan, C. 1999, ApJ, 527, 435 

Tarbell, T. D. 2006, in AAS/Solar Physics Division Meeting, 36.02 

Title, A. M., Tarbell, T. D., Topka, K. P., et al. 1989, ApJ, 336, 475 

Tobias, S. M., Brummell, N. H., Clune, T. L., & Toomre, J. 2001, ApJ, 549, 1183 

Torok, T. & Kliem, B. 2005, ApJ, 630, L97 



25 



Venkatakrishnan, P. 1986, Nature, 322, 156 

Vogler, A. 2003, PhD thesis, University of Gottingen, Germany, 

|http://webdoc.sub.gwdg.de/d iss/2004 /voegler 

Vogler, A., Shelyag, S., Schiissler, M., et al. 2005, A&A, 429, 335 

Weiss, N. O. 1966, Royal Society of London Proceedings Series A, 293, 310 

Zirin, H. 1972, Sol. Phys., 22, 34 

Zwaan, C 1978, Sol. Phys., 60, 213 

Zwaan, C. 1985, Sol. Phys., 100, 397 

Zwaan, C. 1987, ARA&A, 25, 83 



2G 



